% This file take the steady-state the file "../toDynare.mat"
% It writes the dynare file "code_dynare.mod" which is the dynamic model
% to use perturbation techniques
% It doubles check that the steady state is OK, thanks to the "resid" function
% We consider a G shock only


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%        Clearing data         %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
warning('off','all');
isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;
isOctave && warning('off', 'Octave:array-as-logical');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%   Loading steady-state data    %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load ../toDynare % load the file produced by the Julia code


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%    Defining model variables    %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Persistence of the shock
rho_u = 0.99;

% Normalization of the standard deviation for keeping unchanged the NPV of G 
u0 = .01*(1-rho_u/(solution.R))/solution.G; 

% Importing some variables from toDynare.mat and simplifying notation
trunc  = Ramsey.truncatedModel.truncatedAllocation;
lagMul = Ramsey.lagrangeMult;

Ntot  = Ramsey.truncatedModel.Ntot;
Pi_h  = trunc.Pi_h';
S_h   = trunc.S_h;
K     = solution.K(1);
Ltot  = solution.L(1);
alpha = eco.alpha(1);
Y     = K^alpha*Ltot^(1-alpha);
FL    = (1-alpha)*K^alpha*Ltot^(-alpha);
xsis  = Ramsey.truncatedModel.xsis;
ys    = eco.ys;
l     = Ramsey.l;
taut  =  (( 1.0/eco.phi + 1.0)*eco.tau)/( 1.0/eco.phi +1.0 - eco.tau);
x     = trunc.c_h - (1/(eco.chi*(1/eco.phi + 1)))*(trunc.l_h).^(1+1/eco.phi);
omega_h  = Ramsey.weights.omega_h;    %per history
omegab_h = Ramsey.weights.omegab_h;   %per history

% Matrix P of constrained histories
P                 = ones(Ntot,1);
P(trunc.ind_cc_h) = 0;

% Utility function
if eco.sigma==1
    util = @(c) log(c);
else
    util = @(c) (c^(1-eco.sigma)-1)/(1-eco.sigma);
end


% Welfare
iweight = omega_h./S_h;% individual Pareto weights
W = sum(iweight.*S_h.*util(x));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%   Opening the Dynare file    %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

file = 'code_dynare.mod';
fid  = fopen(file, 'w');
formatSpec = '%25.20f';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%     Variables and params     %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['var\n'] ; fprintf(fid, str);
str = ['u r w FK FL A B K mu TFP Ltot Ctot Y W\n'] ; fprintf(fid, str);
str = ['ut Kt Ct Bt wt rt taut tauk mut Zt Yt Lt tkt kappa BoYt G l\n'] ; fprintf(fid, str);

for h = 1:Ntot
    str = ['a',num2str(h),' '] ; fprintf(fid,str);     %  end-of-period
    str = ['at',num2str(h),' '] ; fprintf(fid, str);   % beginning-of-period
    str = ['x',num2str(h),' '] ; fprintf(fid, str);
    str = ['psib',num2str(h),' '] ; fprintf(fid, str);
    str = ['lambdacb',num2str(h),' '] ; fprintf(fid, str);
    str = ['lambdacbt',num2str(h),' '] ; fprintf(fid, str);    % last period
    if mod(h,10)==0; %multiple de 10
        str = ['\n'] ;         fprintf(fid, str);
    end;
end;

str = ['; \n\n'] ; fprintf(fid, str);
str = ['varexo eps; \n\n'] ; fprintf(fid, str);
str = ['parameters\n'] ; fprintf(fid, str);
str = ['beta alpha phi abar delta gamma chi rho_u Gss;\n\n'] ; fprintf(fid, str);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%        Initialization        %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['alpha','   = ',num2str(alpha,formatSpec),';\n']; fprintf(fid, str);
str = ['beta','    = ',num2str(eco.beta,formatSpec),';\n']; fprintf(fid, str);
str = ['phi','     = ',num2str(eco.phi,formatSpec),';\n']; fprintf(fid, str);
str = ['abar','    = ',num2str(eco.a_min,formatSpec),';\n']; fprintf(fid, str);
str = ['delta','   = ',num2str(eco.delta,formatSpec),';\n']; fprintf(fid, str);
str = ['gamma','   = ',num2str(eco.sigma,formatSpec),';\n']; fprintf(fid, str);
str = ['rho_u','   = ',num2str(rho_u,formatSpec),';\n']; fprintf(fid, str); %HERE
str = ['Gss','     = ',num2str(solution.G,formatSpec),';\n']; fprintf(fid, str);
str = ['chi','     = ',num2str(eco.chi,formatSpec),';\n']; fprintf(fid, str);

equa = 1; % Numbering equations
tol = 0;  % threshold for transition, to avoid to consider very small transitions
str = ['\n\n'] ; fprintf(fid, str);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%       Model definition       %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['model;\n\n']; fprintf(fid, str);

%%%%%%%% History-wise FOCs & constraints
for h = 1:Ntot
    
    % Budget constraint
    str = ['x',num2str(h),' = (1/(chi*taut))*l^(1+1/phi)*(',num2str(trunc.y0_h(h),formatSpec),')^taut + (1+r)*','at',num2str(h),'-a', num2str(h), ...
           '; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

    % Definition of atilde (at)
    str = ['at',num2str(h),' = 0 ']; fprintf(fid, str);
    for hi = 1:Ntot
        if (Pi_h(h,hi))>tol   % hi to h
          str = ['+ ',num2str(S_h(hi)*Pi_h(h,hi)/S_h(h),formatSpec),'*a',num2str(hi),'(-1)']; fprintf(fid, str);
        end
    end
    str = ['; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

    % Euler equation
    if P(h)==1 % not constrained
        str = [num2str(xsis.xsiuE(h),formatSpec),'*(x',num2str(h),')^-gamma =', num2str(trunc.resid_E_h(h),formatSpec),' + beta*(1+r(+1))*(']; fprintf(fid, str);
        for hp = 1:Ntot
            if (Pi_h(hp,h)*xsis.xsiuE(hp))>tol
                str = ['+ ',num2str(Pi_h(hp,h)*xsis.xsiuE(hp),formatSpec),'*(x',num2str(hp),'(+1) )^-gamma']; fprintf(fid, str);
            end
        end
        str = ['); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
    else
        str = ['a',num2str(h),' =  ',num2str(trunc.a_end_h(h),formatSpec),';']; fprintf(fid, str);str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
    end
    
    % Definition psi_hat_bar
     str = ['psib',num2str(h),' = - mu*',num2str(S_h(h),formatSpec ) ,' + ',num2str(omegab_h(h)*xsis.xsiuE(h),formatSpec ),'*(x',num2str(h),')^-gamma - (lambdacb',num2str(h),...
            ' - (1+r)*lambdacbt',num2str(h), ')*',num2str(xsis.xsiuE(h),formatSpec ),'*(-gamma*(x',num2str(h),')^(-gamma-1)); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

    % FOC individual savings (at): Euler on psib
    % remark: psib(h) = S(h)*psi(h) = S(h)*sum_hp(psib(hp)/S(hp))
    if P(h)==1
        str = ['psib',num2str(h),'  = ', num2str(S_h(h),formatSpec),'*beta*(1+r(+1))*(0']; fprintf(fid, str);
        for hp = 1:Ntot
            if Pi_h(hp,h)>tol   %be careful transpose from h to hp
                str = ['+ ',num2str(Pi_h(hp,h),formatSpec),'*psib',num2str(hp),'(+1)/',num2str(S_h(hp),formatSpec)]; fprintf(fid, str);
            end
        end
        str = ['); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
    else
        str = ['lambdacb',num2str(h),' =  0; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
    end

    % lambdactc: lambdac previous period
    str = ['lambdacbt',num2str(h),' = 0 ']; fprintf(fid, str);
    for hi = 1:Ntot
        if Pi_h(h,hi)>tol  %be careful transpose
            str = ['+ ',num2str(Pi_h(h,hi),formatSpec),'*lambdacb',num2str(hi),'(-1)']; fprintf(fid, str);
        end
    end
    str = ['; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
  
end
str = ['\n\n %% FOC planner \n\n']; fprintf(fid, str);

% FOC planner l
str = ['0 = ((1 +1/phi)/(chi*taut))*(l^(1/phi))*( ']; fprintf(fid, str);
for h = 1:Ntot
    if S_h(h)>tol   %be careful transpose from h to hp
        str = ['+ psib',num2str(h),'*(',num2str(trunc.y0_h(h),formatSpec),')^taut']; fprintf(fid, str);
        if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
    end
end
str = [') - mu*( ']; fprintf(fid, str);
for h = 1:Ntot
    if S_h(h)>tol   %be careful transpose from h to hp
        str = ['+',num2str(S_h(h),formatSpec),'*(l^(1/phi)/chi)*',num2str(trunc.y0_h(h),formatSpec),'^taut',...
               ' -',num2str(S_h(h),formatSpec),'*FL*',num2str(trunc.y0_h(h),formatSpec),'^(1+taut/(1/phi+1))']; fprintf(fid, str);
        if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
    end
end
str = ['); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% FOC planner B
str = ['mu = beta*(mu(+1)*(1+FK(+1)));  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% FOC planner r
str = ['0 = 0 ']; fprintf(fid, str);
for h = 1:Ntot
    if S_h(h)>0
        str = ['+ at',num2str(h),'*psib',num2str(h),' + lambdacbt',num2str(h),'*',num2str(xsis.xsiuE(h),formatSpec),'*(x',num2str(h),'^-gamma)']; fprintf(fid, str);
    end
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% FOC planner taut
str = ['\n %% FOC planner tau :', num2str(equa),'\n']; fprintf(fid, str);
str = ['0 = ( (l^(1+1/phi))/(chi*taut) )*(']; fprintf(fid, str);
for h = 1:Ntot
    str = ['+ psib',num2str(h),'*(',num2str(trunc.y0_h(h),formatSpec),'^taut)*(-1/taut + log(',num2str(trunc.y0_h(h),formatSpec),'))'];fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end
   str = [' )- mu*(l/(1/phi+1))*( '] ;fprintf(fid, str);
for h = 1:Ntot
    str = ['+',num2str(S_h(h),formatSpec),'*log(',num2str(trunc.y0_h(h),formatSpec),')*(  (l^(1/phi))/(chi)*',num2str(trunc.y0_h(h),formatSpec),'^taut'];fprintf(fid, str);   
    str = [' - FL*',num2str(trunc.y0_h(h),formatSpec),'^(1+taut/(1/phi+1))             )'];fprintf(fid, str);       
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end
str = ['); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Aggregate Constraints

% Public spending shock
str = ['G = Gss*(1+u);\n'] ; fprintf(fid, str); str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Resource contraint
str = ['G + Ctot + K = Y - delta*K(-1) + K(-1);\n'] ; fprintf(fid, str); str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Financial market clearing
str = ['A = B + K; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Production side: factor prices and output. Constant TFP
str = ['FL = (1 - alpha)*TFP*(K(-1)/Ltot)^alpha; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['FK = alpha*TFP*(K(-1)/Ltot)^(alpha-1)-delta; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['u = rho_u*u(-1) + eps;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['TFP = 1; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Y = TFP*K(-1)^alpha*Ltot^(1-alpha);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Aggregate labor supply
str = ['Ltot = l*('] ; fprintf(fid, str);
for h = 1:Ntot
    str = ['+', num2str(S_h(h),formatSpec),'*',num2str(trunc.y0_h(h),formatSpec),'^(1+ taut/(1/phi+1))'] ; fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end;
str = [');   %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Aggregate consumption
str = ['Ctot = 0'] ; fprintf(fid, str);
for h = 1:Ntot
    str = ['+',num2str(S_h(h),formatSpec),'* (x',num2str(h),'+(( l^(1 + 1/phi)) /(chi*(1/phi+1)))*', num2str(trunc.y0_h(h),formatSpec),'^taut )'] ; fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end;
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Wage
str = ['w =  (1/chi)*(1/taut +1/(1/phi+1) ) *l^( (1+1/phi)*(1+1/phi)/(1/phi+1 +taut) )  ;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Taxes
str = ['tauk = 1-r/FK;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['kappa = w/(FL^( (1/phi+1)*taut/(1/phi+1+taut) ));  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Aggregate savings
str = ['A = 0'] ; fprintf(fid, str);
for h = 1:Ntot
    str = ['+', num2str(S_h(h),formatSpec),'*a',num2str(h)] ; fprintf(fid, str);
end;
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

% Welfare
if eco.sigma==1
    str = ['W = 0'] ; fprintf(fid, str);
    for h = 1:Ntot
        str = ['+', num2str(omega_h(h),formatSpec),'*log(x',num2str(h),')'] ; fprintf(fid, str);
        if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
    end
    str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
else
    str = ['W = 0'] ; fprintf(fid, str);
    for h = 1:Ntot
        str = ['+', num2str(omega_h(h),formatSpec),'*( (x',num2str(h),')^(1-gamma)-1)/(1-gamma))'] ; fprintf(fid, str);
        if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
    end
    str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
end;


%%%%%%%% Log-deviation variables
% There are denoted with 't': for variable x, xt = (x -x_ss)/x_ss, where x_ss is the steady state value of x.

str = ['ut = 100*u;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Kt = 100*(K/',num2str(solution.K(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Ct = 100*(Ctot/',num2str(solution.C(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Bt = 100*(B/',num2str(solution.B(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['wt = 100*(w/',num2str(solution.w(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['rt = 100*(r - ',num2str(solution.R(1)-1,formatSpec),');  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['tkt = 100*(tauk -',num2str(eco.tauk,formatSpec),');  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['mut = 100*(mu/',num2str(Ramsey.lagrangeMult.mu,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Zt = -ut;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Lt = 100*(Ltot/',num2str(solution.L,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Yt = 100*(Y/',num2str(solution.Y(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['BoYt = 100*(B/Y*',num2str(solution.Y/((1+eco.tauc)*solution.A-solution.K),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['end;\n\n'] ; fprintf(fid, str);

str = ['%%%%%%%% end of the model part\n\n'] ; fprintf(fid, str);

%%%%%%%% end of the model part


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%     Steady-state model       %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['steady_state_model;\n\n'] ; fprintf(fid, str);
for h = 1:Ntot
    str = ['a',num2str(h),'         = ',num2str(trunc.a_end_h(h),formatSpec),';\n'] ; fprintf(fid, str);
    str = ['at',num2str(h),'        = ',num2str(trunc.a_beg_h(h),formatSpec),';\n'] ; fprintf(fid, str);    
    str = ['x',num2str(h),'         = ',num2str(x(h),formatSpec),';\n'] ; fprintf(fid, str);
    str = ['lambdacb',num2str(h),'  = ',num2str(lagMul.lambdacb(h),formatSpec),';\n'] ; fprintf(fid, str);
    str = ['lambdacbt',num2str(h),' = ',num2str(lagMul.lambdact(h),formatSpec),';\n'] ; fprintf(fid, str);
    str = ['psib',num2str(h),'      = ',num2str(lagMul.psib(h),formatSpec),';\n'] ; fprintf(fid, str);   % recall it is psi multiplied by the size of the bin
end;

str = ['l     = ', num2str(Ramsey.l,formatSpec),';\n'] ; fprintf(fid, str);
str = ['r     = ', num2str(solution.R-1,formatSpec),';\n'] ; fprintf(fid, str);
str = ['w     = ', num2str(solution.w,formatSpec),';\n'] ; fprintf(fid, str);
str = ['FK    = ', num2str(1/eco.beta-1,formatSpec),';\n'] ; fprintf(fid, str);
str = ['FL    = ', num2str(FL,formatSpec),';\n'] ; fprintf(fid, str);
str = ['K     = ', num2str(solution.K(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['Ltot  = ', num2str(solution.L(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['mu    = ', num2str(lagMul.mu,formatSpec),';\n'] ; fprintf(fid, str);
str = ['TFP   = ', num2str(1),';\n'] ; fprintf(fid, str);
str = ['A     = ', num2str(solution.A(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['B     = ', num2str(solution.A(1)-solution.K(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['u     = ', num2str(0),';\n'] ; fprintf(fid, str);
str = ['tauk  = ', num2str(eco.tauk,formatSpec),';\n'] ; fprintf(fid, str);
str = ['taut  = ', num2str(taut,formatSpec),';\n'] ; fprintf(fid, str);
str = ['Ctot  = ', num2str(solution.C,formatSpec),';\n'] ; fprintf(fid, str);
str = ['Y     = ', num2str(solution.Y,formatSpec),';\n'] ; fprintf(fid, str);
str = ['kappa = ', num2str(eco.kappa,formatSpec),';\n'] ; fprintf(fid, str);
str = ['W     = ', num2str(W,formatSpec),';\n'] ; fprintf(fid, str);

str = ['ut    = 0;\n'] ; fprintf(fid, str);
str = ['Kt    = 0;\n'] ; fprintf(fid, str);
str = ['Ct    = 0;\n'] ; fprintf(fid, str);
str = ['Bt    = 0;\n'] ; fprintf(fid, str);
str = ['wt    = 0;\n'] ; fprintf(fid, str);
str = ['rt    = 0;\n'] ; fprintf(fid, str);
str = ['Tt    = 0;\n'] ; fprintf(fid, str);
str = ['tkt   = 0;\n'] ; fprintf(fid, str);
str = ['mut   = 0;\n'] ; fprintf(fid, str);
str = ['Zt    = 0;\n'] ; fprintf(fid, str);
str = ['Lt    = 0;\n'] ; fprintf(fid, str);
str = ['Yt    = 0;\n'] ; fprintf(fid, str);
str = ['BoYt  = 0;\n'] ; fprintf(fid, str);
str = ['G     = Gss;\n'] ; fprintf(fid, str);

str = ['end;\n % end of the steady state part\n'] ; fprintf(fid, str);
%%%%%%%% end of the steady state part


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%   Steady-state computation    %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['resid;\n\n'] ; fprintf(fid, str);
str = ['steady(nocheck);\n\n'] ; fprintf(fid, str);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%   Defining aggregate shock   %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['shocks;\n var eps; stderr ',num2str(u0,formatSpec),';\n end;\n'] ; fprintf(fid, str);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%     Running  stoch_simul     %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['stoch_simul (order=1,irf=500) ut Bt taut kappa tauk Yt Ct Kt Lt mut W r;'] ; fprintf(fid, str);
isOctave && fflush(fid);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%       Launching Dynare       %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

command = ['dynare ',file, ' noclearall'];
eval(command);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%         Saving data          %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

name = ['fig_opti',num2str(rho_u),'.mat'];
save (name,'ut_eps','Bt_eps','taut_eps','kappa_eps','tauk_eps','Yt_eps','Ct_eps','Kt_eps','Lt_eps','mut_eps','W','W_eps')






